Self-Similar Evolution of Cosmological Density Fluctuations 



Bhuvnesh Jain^'^ and Edmund Bertschinger^ 
^Department of Physics, MIT, Cambridge, MA 02139 USA 
^Max-Planck- Institut fiir Astrophysik, Karl-Schwarzschild-Strasse 1, 85748 Garching, 

Germany; bjain@mpa-garclLing.mpg.de 

ABSTRACT 

The gravitational evolution of scale free initial spectra P{k) oc in an 
Einstein-de Sitter universe is widely believed to be self-similar for —3 < n < 4. 
However, for —3 < n < — 1 the existence of self-similar scaling has not 
been adequately demonstrated. Here we investigate the possible breaking of 
self-similar scaling due to the nonlinear contributions of long wave modes. Eor 
n < —1 the nonlinear terms in the Eourier space fluid equations contain terms 
that diverge due to contributions from wavenumber A; — )• (the long wave limit). 
To assess the possible dynamical effects of this divergence the limit of long wave 
contributions is investigated in detail using two different analytical approaches. 

Perturbative contributions to the power spectrum are examined. It is shown 
that for n < —1 there are divergent contributions at all orders. However, at 
every order the leading order divergent terms cancel out exactly. This does not 
rule out the existence of a weaker but nevertheless divergent net contribution. 
The second approach consists of a non-perturbative approximation, developed to 
study the nonlinear effects of long wave mode coupling. A solution for the phase 
shift of the Eourier space density is obtained which is divergent for n < —1. 
A kinematical interpretation of the divergence of the phase shift, related to 
the translational motion induced by the large-scale bulk velocity, is given. 
Our analysis indicates that the amplitude of the density is not affected by the 
divergent terms and should therefore display the standard self-similar scaling. 
Thus both analytical approaches lead to the conclusion that the self-similar 
scaling of physically relevant measures of the growth of density perturbations is 
preserved. 

Subject headings: cosmology: theory — large-scale structure of universe — 
galaxies: clustering — galaxies: formation 
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1. Introduction 

The self-similar scaling of density perturbations with scale free initial conditions in a 
spatially flat universe has been a useful theoretical tool in studying structure formation. It 
has been widely used to study gravitational clustering in cosmology and has been tested 
by several studies using N-body simulations. However self-similar scaling for scale free 
initial spectra P{k) oc k'^ has not been adequately demonstrated for n < —1 because the 
requirements of dynamic range get increasingly difficult to meet as n gets smaller. Indeed 
results of some two dimensional studies suggest a breaking of self-similar scaling for n = —2 
in three dimensions. Analytical analyses have been limited to the observation that the 
linear peculiar velocity field diverges for n < — I, but the linear density contrast does not 
diverge provided n > —3. This would suggest that while there may be formal problems 
with establishing self-similarity for n < — I, in practice it should hold as long as n > —3. 

Our goal in this paper is to analyze the dynamics of the coupling of long wave modes 
by analytical techniques to address the question of whether self-similar scaling is broken for 
scale free spectra with n < —1. A point that will be highlighted in the following sections 
is that the answer can depend on the particular statistic used to pose the question. The 
question of real interest concerns the self-similar growth of measures of density perturbations 
that relate to the formation of structure. Therefore our attempt will be to identify such 
statistical measures and examine their scaling behavior. An analysis of scale free N-body 
simulations that addresses the same issues is presented in a following paper (Bertschinger 
& Jain 1995, Paper II). 

Section I.I gives a detailed discussion of the concept of self-similar scaling and its 
application to cosmology. The statistical nature of self-similar scaling is emphasized and 
an account of previous work that motivated our study is given. This is also done with 
a view to anticipating some of the subtler issues relating to a possible breakdown of 
self-similarity that emerge in our subsequent analysis. Section 2 provides an assessment 
of whether —3 < n < — I is expected to yield self-similar evolution on the basis of simple 
dynamics. We argue that the issue can only be settled by a full consideration of the 
dynamical coupling of long wave modes, rather than by studying the convergence of 
particular statistics using linear solutions. In our analysis we work with the Fourier space 
density field as it quantifies the relative amounts of power on different scales most directly. 
In Section 3 we use perturbation theory to study self-similar scaling. We begin in Section 
3.1 by formulating perturbation theory in a way that obeys this scaling at every order 
provided there are no long wave divergences. In Section 3.2 we demonstrate that there are 
potentially divergent perturbative contributions to the power spectrum, but the leading 
order contributions exactly cancel out. Section 4 presents an alternative, non-perturbative 
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approximation to estimate the coupling of long wave modes. The evolution of the amplitude 
and phase of Fourier modes is studied separately and a solution for the rms phase shift is 
presented. Kinematical effects which do not affect perturbation growth are distinguished 
from dynamical ones which do in order to assess the scaling of the amplitude of the density. 
We conclude in Section 5. 



A physical system is expected to display self-similar evolution if there is no preferred 
scale in the system, either in the initial conditions or in its dynamical behavior. The 
differential equations governing the evolution of such a system generally admit of a 
self-similar solution. Suppose the basic evolution equation is a partial differential equation 
for the phase space density /(x, t), where x is the spatial position, pis the momentum, 
and t is the time variable. In a self-similar system it is possible to re-cast the equation in 
a form with a solution f = t" f{x/t^,p/V), where / is in general an unknown function. If 
the constant coefficients a, /3, and 7 are known then the time dependence of / is present 
only through the rescaled x and p coordinates, aside from the overall factor of t" . This 
special form of the solution is defined to be self-similar: the phase space density at time ^2 
is related to that at time ti as 



where Xi = ^2(^1/^2) , and pi = P2 (^i/^2)"'- Equation (I) explicitly demonstrates that 
the phase space density for any (x^p) at all times ^2 can be obtained merely by re-scaling 
from some chosen time ti. Clearly self-similarity is a powerful constraint because any 
statistical measure constructed from the the phase space density should be described by the 
appropriate scaling of coordinates consistent with equation (I). 

We now consider the similarity properties of gravitational dynamics in a zero-pressure 
Einstein-de Sitter cosmology. An Einstein-de Sitter universe refers to the model with the 
cosmological density parameter = ftmatter = 1 and zero cosmological constant, so that the 
universe is spatially fiat. The gravitational interaction also does not pick a special length 
scale. Eurther let the initial power spectrum be a power law, P{k) oc A;", over length scales 
of interest. Thus so far there is no preferred length scale in the system. The amplitude of 
the power spectrum can be used to define a characteristic physical length scale: the scale at 
which the rms smoothed density contrast equals unity is the conventional choice. To within 
an order of magnitude it is the scale at which over-densities collapse out of the background 
expansion. The presence of this scale does not interefere with self-similar scaling, rather 
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it provides the reference scale required for scaling the spatial variable according to the 
similarity solution. 

The explicit similarity transformation for the single particle phase space density 
t) is described in Section 73 of Peebles (1980). It is also shown that knowing the 
linear solution is sufficient to fix the indices a, /3 and 7 in terms of the spectral index n 
of the initial spectrum. The resulting self-similar scaling of spatial length scales x, and 
wavenumber scales k is: 

oc a(t)2/(3+'^) ; ~ a(t)-2/(3+'^). (2) 

The similarity solution for / is obtained by dimensional analysis of the differential equation 
describing its evolution. Whether or not the solution applies depends on the details of the 
initial conditions. A popular choice for the initial fiuctuations in cosmology is that of a 
Gaussian random field which is statistically homogeneous and isotropic in space. For a 
given realization, the stochasticity of the initial distribution in space precludes the similarity 
solution for / from being valid. Although the statistical properties of the distribution 
are scale free for a power law initial spectrum, the distribution in any one realization is 
not spatially self-similar, or independent of the spatial scale in any meaningful way. This 
means that the phase space density / for a particular realization cannot obey the similarity 
solution. All the results on self-similar scaling that will be discussed below are therefore 
valid only for ensemble averages (averages across different independent realizations). This 
distinguishes the self-similar scaling of cosmological perturbations from standard examples 
of self-similar systems. 

The ensemble averages of / or products of / do evolve self- similarly because ensemble 
averaging removes the stochastic character of the initial conditions. Thus the self-similar 
scaling for the 2— point correlation function ^{x^t) can be obtained from the formal solution 
for /, and it is valid even though the self-similar solution for / does not hold. Its validity 
can be verified by dimensional analysis of the evolution equation for ^ that is obtained by 
taking moments of the BBGKY hierarchy equations (e.g., equation 1.1 of Peebles 1980). 
The solution for the power spectrum is obtained by Fourier transforming (^(a;,t), and is 

P{k,t) = a'^"k-'^P{ka"/ko) , (3) 

where a = 2/(3 -|- n), ko is a constant which must be determined from the initial conditions, 
and P is an unspecified dimensionless function. It is easy to verify that the linear spectrum 
Pii(A;,t) oc a^k^ is consistent with this functional form. 

Likewise the scaling of all statistical measures defined as ensemble averages of products 
of / (and of their momentum moments) can be straightforwardly determined. Using the 
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ergodic theorem the solutions for ensemble averages can then be applied to averages over 
sufficiently large volumes in space. A spatial statistic which follows the self-similar solution 
is a function of the spatial variable scaled by a power of time, rather than of time and 
space separately. This provides for a self-similarity in time (in this statistical sense) in the 
evolution of structure. 

The discrete nature of particles (N-body or galaxies) introduces a scale, namely the 
mean interparticle separation, which breaks the idealized self-similar scaling of a perfect 
fluid. Such a departure from perfect self-similarity is typical of all realistic physical 
systems. The notion of intermediate asymptotic self-similarity, i.e., self-similar scaling 
over a restricted range of parameters, is used in such situations (Barenblatt 1979). In the 
cosmological context it simply means that the range of length scales over which self-similar 
scaling is accurately followed are restricted to be sufficiently larger than the interparticle 
separation (and in the case of N-body simulations, the force resolution scale). 

Intermediate asymptotic self-similarity is a useful property even for realistic 
cosmological spectra like the CDM spectrum which are not scale free. The physical 
processes at work in the radiation dominated era imprint characteristic length scales on the 
spectrum. These spectra are nevertheless approximate power laws on a restricted range 
of A;, over which their evolution may be well described by the similarity solution for the 
corresponding scale free spectrum. Thus the CDM spectrum has n ~ —2 on galactic scales 
and n ~ —1 on cluster/supercluster scales; therefore the study of scale free spectra with 
— l^n^ — 2 is relevant for understanding the development of large scale structure in 
a CDM-like model. Moreover, recently the idea of self-similar scaling has been extended 
to provide a more direct scaling description of the evolution of CDM-like initial spectra 
(Hamilton et al. 1991; Peacock & Dodds 1994; Mo, Jain & White 1995). 

An aspect of self-similarity which merits attention is the range of n, the spectral index 
of the initial spectrum, for which the statistics characterizing the growth of perturbations 
are well defined. More precisely, n must be restricted from below and above to prevent 
statistical measures of interest from diverging as the size of the system is made infinitely 
large and the interparticle spacing made infinitesimally small, respectively. In a finite 
system such a divergence is manifested by the infiuence of the largest (or smallest) scales on 
the evolution of all intermediate scales of the system. This occurs if either the statistic used 
is ill defined even in the initial configuration, or the dynamical infiuence of increasingly 
large or small scales is unbounded. If the former is true in an otherwise reasonable initial 
configuration, then it must mean that the particular statistic is not a suitable measure of 
the properties of the system (similar to the case of well behaved probability distributions 
having divergent moments). However, if the breaking of self-similarity is due to a divergent 
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dynamical effect in a statistic of interest tlien it bears closer examination. The goal of 
this work is to examine the possible breaking of self-similar evolution for power law initial 
spectra with a view to assessing its influence on the formation of structure. 

Early studies of self-similar evolution in cosmology include those of Peebles (1974); 
Press & Schechter (1974); Davis & Peebles (1977); and Efstathiou & Eastwood (1981). 
Davis & Peebles (1977) made a detailed analysis of the BBGKY hierarchy equations 
and presented approximate solutions for the deeply nonlinear regime based on the stable 
clustering ansatz. Efstathiou et al. (1988) tested self-similar scaling in N-body simulations 
of scale free spectra with n = —2, —1,0, 1. They examined the scaling of the correlation 
function (^(a;,t), and of the multiplicity function describing the distribution of bound 
objects. They verified the predicted scaling for both statistics, and found consistency with 
the picture of hierarchical formation of nonlinear structure on increasingly large length 
scales. Their results for n = —2 did not match with the self-similar scaling as well as the 
other cases. Bertschinger & Gelb (1991) used better resolution simulations to address these 
questions and also found similar results. These authors concluded that the reason for the 
weakness of the n = —2 results was the finite size of their simulation box, as the n = —2 
case has more power on large scales and therefore requires a larger box-size to approximate 
the infinite volume limit with the same accuracy as larger values of n. 

Recently Lacey & Cole (1994) have examined the self-similar scaling of the number 
density of nonlinear clumps for scale-free spectra. Their results indicate that self-similar 
scaling works reasonably well for the statistics they measure, even for n = —2. Mo, Jain & 
White (1995) reach the same conclusion for the correlation function and power spectrum. 
However, in all the above studies the self-similar scaling of an averaged quantity is tested 
without making a comparison with alternate scalings. This makes it difficult to distinguish 
the effect of limited numerical resolution from a real breakdown of self-similar scaling. 
In Paper II we demonstrate that even with the current state-of-the-art simulations it is 
extremely difficult to use an averaged quantity like the power spectrum to critically test 
self-similar scaling for n = —2. 

The N-body results of Ryden & Gramann (1991), and Gramann (1992) suggest that 
the n = —2 case is different for a fundamental reason. They studied n = —1 simulations 
in two dimensions, which are the analog of n = — 2 in three dimensions, and examined 
the scaling of the phase (Ryden & Gramann 1991), and then both phase and amplitude 
(Gramann 1992) of the Eourier transform of the density field. The scaling was found to be 
different from the standard self-similar scaling. Characteristic wavenumber scales, instead of 
following the self-similar scaling, given in two dimensions by kgsit) oc a(t)~^/(^"'"") oc a(t)~^, 
showed the scaling k oc a{t)~^ . Other studies in two dimensions also suggest that a 
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transition in nonlinear evolution occurs at n = — 1 (Klypin & Melott 1992, and references 
therein). Motivated by Gramann's results we had re-examined the n = —2 simulation 
presented in Bertschinger & Gelb (1991) and found that the results were ambiguous, and 
that a bigger simulation would be needed to provide a definitive answer. This had provided 
the initial motivation for our analytical investigation as well. 

It has been noted all along in the literature that the formal bounds on n for the 
self-similar solution to be applicable are — 1 < n < 1. The requirements of an upper (lower) 
limit are made to prevent the single particle velocity dispersion from diverging due to 
contributions from small (large) length scales. These bounds on n are clearly stated as the 
domain of applicability of self-similar scaling in Peebles & Davis (1977), Efstathiou et al. 
(1988), and in the recent review of Efstathiou (1990). However, it appears to be implicitly 
believed that self-similar scaling is applicable for n > —3, rather than n > —1. This is 
probably because the divergence of the bulk velocity field need not affect the growth of 
perturbations. The primary quantity that measures perturbation growth is the rms density 
contrast which is indeed convergent for n > —3 as A; — )• 0. Thus Peebles (1993, p. 545) 
presents the standard self-similar scaling as being applicable for —3 < n < 4 (increasing 
the upper limit to n = 4 relies on the asymptotic behavior of second order contributions 
to the density). Efstathiou (1990) was more cautious, but nevertheless hoped that: "If n 
lies outside this range (i.e., —1 < n < 1), the clustering may still approximate self-similar 
evolution over restricted ranges of length and time, although n > —3 is required to ensure 
that clustering proceeds from small to large scales." 



We suppose for simplicity that the matter distribution after recombination may be 
approximated as a pressureless fiuid with no vorticity. We further assume that peculiar 
velocities are nonrelativistic and that the wavelengths of interest are much smaller than 
the Hubble distance cH~^ so that a nonrelativistic Newtonian treatment is valid. Using 
comoving coordinates x and conformal time dr = dtj ait), where ait) is the expansion scale 
factor, the nonrelativistic cosmological fiuid equations are 



2. Long Wave Divergences for n < —1 
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where a = da /dr. Note that v = dx/dr is the proper pecuhar velocity, which we take to 
be a potential field so that v is fully specified by its divergence: 

e = V-v. (5) 

We assume an Einstein-de Sitter (0 = 1) universe, with a oc t'^^^ oc r^. We will also assume 
that the initial (linear) density fiuctuation field is a gaussian random field. 

To quantify the amplitude of fiuctuations of various scales it is preferable to work with 
the Fourier transform of the density fiuctuation field, which we define as 

HKr) = J^^e-^'-^Six,r), (6) 

and similarly for 6{k,T). The power spectrum (power spectral density) of 6[x^t) is defined 
by the ensemble average two-point function, 

{6{h,r)6{k2,r)) = P{h,T)6G{h + k2) , (7) 

where Sjy is the Dirac delta function, required for a spatially homogeneous random density 
field. For a homogeneous and isotropic random field the power spectrum depends only on 
the magnitude of the wavevector. The contribution to the variance of 6[x^ r) from waves in 
the wavevector volume element d^k is P[k^T)d^k. 

Fourier transforming equations (4) gives: 
38 ^ f k ■ ki 



— + = -J d'kr e{kr,T) 6{k -kr,T) , (8a) 

^ + -e + %8 = - ( d''krk^^^^^i^^e{h,T)e{k -h,T) . (8b) 

Ot a J 2kj\k - ki\^ 

The fields on the left-hand side are all functions of k and r. The nonlinear terms on 
the right-hand side of the above equations represent the coupling of modes at all pairs of 
wavevectors (A;i, k — ki)^ which infiuence the evolution of 6 and 6 at the fixed external 
wavevector k. 

In order to study the limiting behavior as one of the pair of wavevectors {ki^k — ki) 
approaches in magnitude (i.e., as the wavelength A = 2Tr/k is made infinitely large), 
consider the variance of the nonlinear terms. For simplicity we take A;i — )• in the integral 
on the right-hand side of equation (8a). Approximating [k — ki) by A;, dropping the 
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dependence on r, and denoting the resulting variance by V'(^) we obtain: 



,l;{k) ~ { ld'kr—^e{kr)8{k) d' k^ —j^ 9* {k^) 8* {k) 



d'kr j d'k2 -—^ (eih) r (4) sik)s*ik) ) . (9) 



k ■ ki k ■ k2 

k-^ k 

Now we make the further approximation of taking ki small enough that the linear solutions 
are valid, thus giving: 

6{h,r)c:,a{T)6^{h) ; e{h,r) c:, -h{T)6^{h) ■ (10) 

Substituting the expression for 6{ki,T) in equation (10) into equation (9) and evaluating 
the ensemble average using the properties of Gaussian random fields we finally obtain: 

i^{k)c^h^8^{<d)P{k) jd^'kri^^^^^ Prr{kr)c^h^8^{<d)P{k)k^ jdkrPrrikr). (11) 

Note that in expressing the 4— point moment of equation (9) in terms of P{k) and P\\{ki) 
we have assumed that 8{k) is a Gaussian random field as well. It is straightfoward to 
demonstrate that the right-hand side of equation (8b) takes the limiting form shown in 
equation (11) as well. 

Equation (11) indicates that if P\\{ki) oc k'1 with n < —1, then the right-hand side 
diverges due to contributions from low k^. Thus a simple examination of the nonlinear 
terms in the cosmological fiuid equations by substituting the initial distribution of the 
density and velocity field demonstrates the possibility of long wave divergences. These 
divergences can potentially be present in solutions for 8{k,T) and 6{k,T) obtained from 
these equations. It is not possible to make any definitive statements, however, because these 
are two coupled differential equations — it is necessary to first separate the equations for 8 
and 6, and then identify the nonlinear terms that affect the amplitude and phase of each 
quantity (since they are complex variables) to determine whether the divergent terms affect 
a particular statistic of interest. This is done in two different ways in Sections 3 and 4. 

Before proceeding with a formal analysis of the divergent nonlinear pieces, we make the 
connection between the divergent nonlinear terms to the advective [v ■ V) terms in the real 
space equations. By tracing back the nonlinear terms on the right-hand side of equations 
(8) to the fiuid equations in real space it can be seen that the terms which contribute to 
equation (9) arise from the v ■ V(5 and v ■ terms. It is easy to see why such terms should 
diverge by examining the relation of the power spectrum of the peculiar velocity to that of 
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the density in hnear theory. Using the hnear solutions of equation (10) and the definition 
6[x^ r) = V • u(x, r), gives, 

Pu.{k,T) = h'Pu{k)/k\ (12) 

where Piiv{k) is the linear power spectrum of the peculiar velocity. The rms bulk velocity 
on a scale x ~ k~^ , Vb{x,T), is given by integrating Pv{k) over k with a window function 
W{kx) (just as one integrates over P{k) for the rms smoothed density contrast): 

Vb{x,Tf = j d^kPii^{k,T)W^{kx) = a^TT j dkPii{k)W^{kx). (13) 

Since W{kx) — )• 1 as A; — )• (see for example the standard top-hat window function), the 
integral on the right-hand side of equation (13) diverges at low k for n < —1 in the same 
manner as the integral in equation (11). Thus via the advective [v ■ V) terms in the fiuid 
equations, the divergence of the nonlinear terms demonstrated in equation (11) can be 
traced to the bulk velocity field on a given scale receiving divergent contributions from 
A; — )• 0, i.e., from the long wavelength modes. 

We can now understand why this divergence may not affect self-similar scaling: the 
bulk velocity field does not in general have any infiuence on the growth of perturbations 
on small scales. In particular, large contributions to the bulk velocity field from long wave 
modes correspond to an almost uniform translation of the fiuid, and therefore should not 
couple to the evolution of 8 at all. This reasoning, and indeed the entire analysis of this 
section, relies on making plausible connections of linearized statistics for 8 and v to their 
nonlinear dynamics. Therefore, while it provides a useful guide to one's intuition, it does 
not substitute for a rigorous examination of the nonlinear dynamics. 

3. Self-Similarity and Perturbation Theory 

3.1. Formalism 

The basic formalism of perturbation theory used here most closely follows that of 
Goroff et al. (1986); it is described in detail in Jain and Bertschinger (1994). We begin by 
writing the solution to equations (8) as a perturbation series, 

oo oo 

8{k,T) = Y,a\T)8i{k) , e{k,T) = Y,h{T)a'-\T)ei{k) (14) 

It is easy to verify that for 1 = 1 the time dependent part of the solution correctly gives 
the linear growing modes 8i oc ^(t) and 6i oc a and that the time-dependence is consistent 
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with equations (8) for all /. To obtain formal solutions for the k dependence at all orders 
we proceed as follows. 



Substituting equation (14) into equations (8) yields, for / > 1, 

ISiik) + 9iik) = Mk ) , 3Siik) + (1 + 2l)9iik) = Bi{k) 

where 

Ai{k) = - 



Bi{k ) = - d'kr / d'k2 8j,{kr + k^ - k) 



kl 



m = l 



l-l 



^ 0,m{kl) Ol-,rn{k2) 



2 m = l 



Solving equations (15) for 8i and 9i gives, for / > 1, 



^^^^^^_ {l + 2l)Mk)-Bi{k) ^ -3A{k) + lB,{k) 



(2/ + 3)(/-l) 



(2/ + 3)(/-l) 



(15) 

(16a) 
(16b) 

(17) 



Equations (16) and (17) give recursion relations for 6i[k) and 0i{k)^ with starting 
values 6i[k) and 6i = —6i. The general solution may be written 



Si[k ) = J d qi - ■ ■ J d qi 6n{qi H V qi - k) F/(^i, ...,qi) 8i{qi) ■ ■ ■ 8i{qi) , 

9i(k) = - j d^qi ■■■ j d^qi6r,{qi H \- qi -k) Gi{qi, . . .,qi) Si{qi) ■ ■ ■ 6i{qi) . 

From equations (16)-(18) we obtain recursion relations for Fi and G;: 

' 7 ^m) 



\ (2/ + 3)(/- 1) 



k ■ kl 

(1 + 21) Ff_m(gm+i, 



Gl-m{(lm + l, ■ ■ ■ T^ll) 



(18a) 



(18b) 



(19a) 



Gi{qi, ...,qi) 



\ ^ Cm ( *?! 7 • • • 7 ([m ) 

h (2/ + 3)(/-l) 
k\k ■ K 



k ' k\ 

3 1^ -r)_m(<?m + l7 • • • 7 



qi) 



+ 1 



Ivy 1^2 



-Gi-miq, 



m + 1 7 



*) 



where ki = qi + ■ ■ ■ + q^, k2 = qm+i + ■ ■ ■ + qi, k = ki + k2 and Fi = Gi = 1. 



(19b) 



To calculate the power spectrum we shall prefer to use symmetrized forms of Fi and G;, 
denoted f/*-* and gJ*'' and obtained by summing the /! permutations of Fi and Gi over their 
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/ arguments and dividing by /!. Since the arguments are dummy variables of integration 
the symmetrized functions can be used in equations (18) without changing the result. The 
recursion relations in equations (19) may be used to compute the power spectrum at any 
order in perturbation theory. Substituting equation (14) into equation (8), we have 



P{k,T)8j,{k + k') 



{8{k,T)8{k\T)) 

a2(T)(,5i(^) 8S')) + o\ 



{8^{k)8s{k')) + {82{k) 82{k')) 



+ {8,{k)8,{k')) 



+ 0(8^: 



l-l 



m = l 



(20) 



1=2 



with / restricted to being an even integer. Equation (20) explicitly shows all the terms 
contributing to the power spectrum at fourth order in the initial density field 8i. To 
distinguish the different terms that contribute to the power spectrum at a given order we 
have introduced Pm,/-m(^), the mth contributing term at order /, defined as: 



{8^{k) 8i_^{k')) = P^^i_^{k) 8Y^{k + k') . 
Substituting for 8i and 8ra in equation (21) gives: 



f211 



Pm,i-mik, t) = a' J d^qi . . . d^q(i_2) ('^1(91) . . . 8i{k - q^ - . . . - q^-i) 8i{q^) . . . 

8i{k - qm - . . . - q{i-2))) Mi(^i, . . . , q(i-2), k) , (22) 

where Mi is a dimensionless function of F^"^ and Fi^^. 

We can now proceed to demonstrate that in the absence of long wave divergences 
(i.e., if the contributions to P{k) from low q are convergent), a perturbative expansion 
can be consistently defined such that the self-similar scaling of equation (2) is obeyed. At 
sufficiently small scales 8pl p > 1 even at the earliest times — hence the perturbative 
expansion breaks down at these scales. This means that even in the absence of a divergence 
as k ^ 00, a high-A; cutoff, A;„, must be used to truncate the perturbative integrals. The 
requirement of a cutoff restricts the nonlinear effects that can be studied perturbatively. 
Nonlinear contributions from wavenumbers q > ku cannot be evaluated, but the contribution 
to any k (from all q < ku) are calculable. We define the high-A; cutoff to be time dependent 
such that: ku{T) oc ksgi^T) oc a"^/^"^"'""^ Once ku{T) is chosen to scale with ksgi^T)^ there is no 
other scaling in the problem, hence it should come as no surprise that the power spectrum 
takes the self-similar form of equation (3). 
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Without loss of generality, consider the contribution to P{k,T) from the term 
Pra,i-m{k,T) defined in equation (21). On taking the ensemble average, the (/ — 2) 
independent phase factors contained in the functions 6\{qi) must cancel pairwise for the 
result to be non-zero (recall that the (5i's are taken to be independent Gaussian random 
variables). Thus we obtain (/ — 2)/2 Dirac-delta functions which reduce the number of 
integration variables to (/ — 2)/2. There are also //2 powers of Pii[qi) = Aq^ present. 
Collecting the relevant factors which provide the k and r dependence, and imposing the 
high-A; cutoff ku{T), we obtain 

P^,i-^{k,T) = a'k^^'-'^l'+-'l'A'l'M,{k^{T)lh), (23) 

where M2 is another dimensionless function. Taking ku{T) = koa~'^^^^~^'"'\ where ko is a 
constant, and introducing a new dimensionless function Ms, we finally obtain 

Pm,l-m{k,T) = a6/(3+-^)^-3^^(^^2/(3+.)/^^^_ (24) 

Equation (24) gives the desired self-similar form for the power spectrum (equation 3), with 
characteristic scales kr (r) oc a-2/(3+'^) in a greement with equation (2). 

Thus at every order in perturbation theory the self-similar scaling of the power 
spectrum, and therefore of physical measures of perturbation growth such as the smoothed 
density contrast in real space constructed from it, is preserved. However, this scaling is 
broken if the perturbative integrals diverge as A; — )• 0, thus requiring a low- A; cutoff. This 
possibility is considered next. 

3.2. Long Wave Divergences in Perturbative Contributions 

The perturbative formalism presented above was used in an earlier work (Jain & 
Bertschinger 1994) to study the second order power spectrum. It was demonstrated that 
at second order in the power spectrum there are terms that are divergent for n < —1 
due to the contribution from A; — )• 0. However, as first shown by Vishniac (1983) the two 
contributing terms, P22 and P13, exactly cancel each other at leading order as A; — )• thus 
keeping the net contribution finite for n > — 3. This cancellation does not prove there is 
no divergence in the power spectrum. We must investigate higher-order terms 6i{k ). It is 
tedious to evaluate the full expressions for 81 for / > 2 and then to form the power spectrum 
contributions Pm,i-m{k). However, we do not need the exact nonlinear power spectrum 
if we are interested only in determining whether leading-order long wave divergences are 
canceled. In this case, it is sufficient to work from the outset with only the leading-order 
divergent parts of 6i[k ). 
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Iterating equations (16) and (17), one finds that the leading-order divergences arise 
from the term with m = 1 in equations (16), with the contribution doubled in equation 
(16b) because of the term m = I — 1. The leading-order divergent contributions are then 

Mk)^8i_S)ak) . Bi{k)^ei_S)ak) . ak) = j d\^8r{q) . (25) 

The leading-order divergence appears at (j' = in the function ({k). Using equation (17) 
and iterating gives the leading-order divergences of 6i and 6i: 

S!ik)r.e,ik)r.tlitls,(k). (26) 



From equations (21) and (26) we arrive at the leading-order divergent part of 



ra.i — ra • 



P-r.,l-ra{k ) ~ . . ^ ^ Pll{k) , (27) 

[m — L)l [I — m — Ij! 



where 

(/-2)/2 



28) 



3q 

The net contribution to the leading-order divergent part of the nonlinear power spectrum 
(20) is 

Now, the sum over m is just the binomial expansion of (1 — 1)'~^. Therefore, the sum 
vanishes for n > 2 and the leading-order divergences cancel at every nonlinear order of 
perturbation theory! 

This surprising result does not prove that P{k^ r) is finite, however. Equation (26) gives 
only the most divergent contribution to Si{k ), ('~^. Terms diverging as ('~'^ or more slowly 
have been neglected. The nonlinear power spectrum may still have an /-point contribution 
that diverges as (C'~^)- The lowest order at which such sub-dominant divergences would 
appear is / = 6: P = a^[Pi5 -\- P24 + -P33 + P42 + Psi)- Using equations (18a) and (21) gives 
the integral expressions for the contributing terms: 

Prs{k) = l5Prr{k)Jd^q,P,,{q,)Jd^q2Prr{q2)F^'\q,,-q,,q2,-q2,k), (30) 

P2,{k) = l2 j d^q, Pn(9i) / ^2^11(92) Pii(|^ - 92!) Fi'\q,, -q,, q^, k - q^) 

xF^'\-q2,q2-k), (31) 
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P33{k) = 9Pii{k) I cfqr Prriqr) j cfq^ ^11(92) Fi'\qr, -q^, k) Fi'\q2, -q^, -k) 

+6 j d^qiPiiiqi) J d^q2 ^11(92) Piii\k - qi - 92!) 

xF3^*''(^i, ^2, k - qi - ^2) F^'\-qi, -^2, qi + q2 - k) . (32) 

The factors of 15, 12, 9 and 6 in equations (30)-(32) come from the number of equivalent 
graphs obtained by relabehng the internal wavevectors, assuming that F(^) and G'(^) are 
fully symmetric in all their arguments. By examining the form of the sub-dominant 
divergent parts of the contributing terms, it appears that the second term in P33 must 
cancel with the divergent part of P24 as (j'l — )• 0, and the other terms shown must cancel 
separately if there is to be a net cancellation. With increasing /, the full expression for 
F;((fi, . . . , (f;) rapidly becomes unwieldy. We did not complete this calculation owing to the 
computational complexity involved, and the fact that a null result would leave open the 
possibility of sub-dominant divergences at still higher orders. 

The results from the analysis of perturbation theory are therefore not conclusive. The 
cancellation of leading divergences is certainly suggestive of an underlying kinematical 
effect which appears in the power counting assessment of the divergence, but cancels out 
on computing the net dynamical influence on the power spectrum. We will interpret this 
cancellation by examining the phase of 6[k) in Section 4. However it is not feasible to 
evaluate all the divergent terms at arbitrary order in perturbation theory, therefore we 
pursue a somewhat different approximation to evaluate long wave mode coupling in the 
next section. 

4. Analytic Approximation for Long Wave Mode Coupling 

The approach in this section relies on assuming that the nonlinear terms in the 
Fourier space cosmological fluid equations (8) are dominated by the coupling of long wave 
modes. With this ansatz the mode coupling contribution is estimated and then checked 
for self-consistency. This allows us to obtain a leading order solution for the phase shift 
as described in Section 4.1. To make further progress we need to make the additional 
assumption that at low A;, the Fourier space density and velocity fields are continuous and 
therefore amenable to a Taylor series expansion. This analysis is presented in Section 4.2, 
and its limitations are discussed. 



4.1. Solution for the Phase Shift 
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In equations (8) the integrands on the right-hand side involve products of 6 and 6 
evaluated at ki and [k — ki). Let 

S{k - ki) = 6{k) + e{k, ki) ; e{k - ki) = e{k) + uj{k, h) , (33) 

where e and io are unknown functions. In this section we shall use "function" to refer to 
random valued fields as well. We shall also suppress the time dependence of 8{k,T) and of 
9{k, t) for convenience, though when we introduce the linear solutions the r dependent part 
will be explicitly written. Substituting equation (33) into equations (8) gives, 

^^^^^ + e{h) = - ( d^kr e{h) \6{k) + e{k, h)] = A{k), (34a) 



^^(^) + 1 0(k) + 1 S{k) = - [ d'k^ k' ^' ' ~ \e{k) + io{k, h)] = B{k). (34b) 



dT a ' ' ' ' J 2kl\k-ki 

In order to estimate the nonlinear effects of long wave modes we assume that the 
integrands on the right-hand side of equations (34) are dominated by the contribution 
from ki <^ k. We then approximate 9{ki) by the linear solutions given in equation (10): 
6i[ki^T) = —a6i[ki)^ because for ki <^ k the amplitude of the density perturbations is 
taken to be very small. Thus we write the right-hand side of equations (34) as 

A[k) = a6[k) Jd'^ki—^ — '^i(^i) +a J'^'^^i~~^ — '^i(^i) e(^, ^i), (35a) 

B{k) = h 0{k) j d^k^ +dj d^k^ ^^8^{k^)u{k, ki). (35b) 

In the expression for B[k) we have multiplied the right-hand side by 2 to include the 
contribution from [k — ki) ^ as required by the symmetry of the integrand. We have also 
explicitly written out the r dependence of 0{ki)^ so that 6i[ki) does not depend on r. We 
now define the integrals: 



a = —i I d^ki -p^8i{ki) , (36) 

where i = \/— T; and. 



E{k) = h I d^ki^^Si{ki)e{k,ki) ; Wik) = h j d^k^^^6^{k^)Lo{k,k^) . (37) 



^-^sSMkM) ■ W{k) = h jd^kj^^ 

Using these definitions equations (34) can be written as 

d6{h) 



dT 



+ e{k) = idk-a6{k) +E{k), (38) 
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'^^ + -e{k) + -^6{k) = ihk-ae{k) +W{k). (39) 




The above equations are exact aside from using the hnear solutions for 6{ki) in 
the right-hand side of equations (34). We have defined a in equation (36) so that it is 
purely real. This can be verified by using the relation of 6i to its complex conjugate: 
'^i(^i) = ^i{~ki)^ which is required to ensure that 6[x) is real. The reason for introducing a 
is that it is independent of k and r, therefore, given the initial density it can be treated 
as a numerical constant. 

We now turn to the issue of long wave divergences. The variance of a is: 

{a') = {a -a) = J d^k^^^^ = iir J dk^Puiki). (40) 

Thus (a^) is a divergent integral for n < —1. To proceed further we need to estimate the 
degree of divergence of the integrals E[k) and W{k). We do so by using equation (33) to 
substitute for e and io in E and W. The resulting expression for the variance of E is: 

mk)\') = Jd'k,Jd'k,(^^ 

x{S^{h)Sl{h) [S{k-h)-S{k)] [s*{k-h)-S*{k)]). (41) 

To simplify this expression we assume that for the purpose of assessing the degree of 
divergence, all the fields involved are well approximated by the linear solution. Then the 
expectation values can be evaluated using the properties of Gaussian random fields. Of the 
twelve terms that result, the leading contribution in the long wave limit arises from the 
term with {Si[ki)6l[k2)) {S[k — ki)6*[k — k2)) in the integrand. This contribution is: 

{\E{k)\^) ^^-h" P{k)8D{^)ea\ (42) 

The variance of the first term on the right-hand side of equation (38), iaa ■ k6[k) is exactly 
the same as the above result for (|£'(A;)p), hence both terms must be retained at the same 
order in evaluating the long wave contribution. Likewise, it is easy to show that W in 
equation (39) is of the same order as the first term on the right-hand side, and is also 
proportional to a in its degree of divergence. 

Equations (38) and (39) can be re-written as a pair of second order differential 
equations in r for 6 and 6. For 6 the result is, with a = k ■ 



S + S { -2ih & + -^ +S (^-h^ &^ - 3id& - -^^ - -E + ih&E + W - E = 0. (43) 
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Since, to leading order, we know the degree of divergence of the variances of all the terms 
involved in this equation, we are now in a position to evaluate the effect of long wave 
divergences. The variables in equation (43) are complex, hence it can be simplified further 
by separating the real and imaginary parts. To this end we express 6 in terms of its 
amplitude and phase as: 

6{k,T) = A{k,T) e"^^^'^K (44) 

For convenience, we further define E = Ee~^'^ , W = We~^'^ , and E = Ee~^'^ . With these 
substitutions equation (43) separates into its real and imaginary parts (respectively) as: 



A + -A + A + 2^aa - a^a^ - ) + Re 



-E + ih&E + W - E 
a 



A [2ct> - 2aa) + A ( ^ + - 3aa ) + Im 



-E + ih&E + W -E 
a 



, (45a) 



. (45b) 



"Re" and "Im" denote the real and imaginary parts, respectively, of the expressions in the 
square brackets. 

We now make the ansatz that (f> ~ 0(a), and that A ~ (9(a''), where < p < 1. Since 
E W ^ 0(a) (in an rms sense) from equation (42), keeping terms of 0{a^) in equation 
(45a) gives, 

A (-(^2 + 2<j)aa - a^c?^ + Re [laaE - e\ = . (46) 

As we shall see below, retaining the term E is required for consistency. We make the 
assumption that at leading order in a the two parts of equation (46) in brackets vanish 
separately (this will also be justified below). The first part gives a quadratic equation for ^, 

- 2^aa + a^c? = , (47) 

which has the solution, (f> = aa = ak ■ a. Thus the leading order solution for (f> is: 

(P{k,T) = a{T)k ■ a + (p,{k), (48) 

where (f>i{k) is the value of the phase at the initial time. 

The solution of equation (48) can be used to justify the assumptions that have been 
made. Firstly, (f> is indeed of 0(a), as assumed at the outset. Further, equations (37) and 
(48) can be used to simplify the expression for E and thereby justify setting the first part 
of equation (46) to separately. To start with let us write E in terms of A and (f>: 

E = h j d^k^^^A^{k^)e'*^^^'^ [A{k - h)e"^^^-^'^ - A{k)e"^^^^] . (49) 
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Differentiating equation (49) with respect to r and multiplying by e ''^^'^'^ gives, 

E = -^ + a /j3A;ii-|iAi(^i)e''^(*^i) [|A(^- ^i) + zA(^- ^i)(^(^- ^i)| 
a J k( 

^^,4,(k-h)-^Hk) _ 1^^^^ ^ iA{k)^{k)}]. (50) 

The leading order terms on the right-hand side above are the two terms with (f>: they are 
at least of 0{a^). However, by substituting <f> = aa into equation (50), and comparing 
with the expression for E that follows from equation (49), it can be seen that these leading 
order terms exactly cancel the contribution from iaaE in equation (46). Thus the surviving 
terms in the second part of equation (46) are all of lower order than 0{a^) — therefore 
they can be neglected in comparison to the first part of the equation which was used to get 
the solution for (f> of equation (48). This establishes the consistency of the approximations 
used to obtain this solution. 

The variance of the phase shift given by equation (48) is: 

airfk' jd'k^^^^P,,{k,), (51) 

where 8cf){k,T) = cf){k,T) — (f>i{k)^ and k and ki are unit vectors. Since A; is a fixed external 
vector the angular integral can be performed so that the final result depends only on the 
magnitude of k: 

= ^airyk' Jdk,P,,ik,). (52) 

Thus the leading order solution for 6[k) involves a growing (and, for n < — 1, divergent) 
phase shift, but there are no contributions to the amplitude at this order. The above 
analysis can be repeated for the velocity divergence 0{k) to verify that the leading order 
result for 0{k) is the same, with equation (48) giving the solution for its phase as well. 
These results were obtained by retaining terms of 0{a^). Since divergent terms of 0(a) 
are also present in the equations we cannot say anything conclusive about the degree of 
(possible) divergence of the amplitude A[k). In the following section we shall address this 
question by expanding the equations to next order in a with some additional assumptions. 

The solution (48) for the phase shift has a simple physical interpretation. As noted in 
Section 2, the linear bulk velocity vj, diverges due to contributions from long wave modes 
(equation 13). The limiting form of the integral given in equation (13) for u^, and that 
of equation (52), is the same. The connection between them can be made more precise 
by imagining a single sine-wave density perturbation in real space: 6[x^t) = 6osin[k ■ x). 
Now suppose that the fiuid in which this perturbation is made is moved with a uniform 
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translational velocity of magnitude U6(t) given by (13) (the scale x in equation 13 has no 
connection to the spatial variable x used here). The distance moved by each fluid element is 
/ dT V}j{t) = V}j{t) a{T) I a{T) , where we have used V}j{t) oc a(T). If the coordinate frame is 
kept fixed relative to this translational motion, then the density perturbation will acquire the 
following time dependence due to the bulk velocity: 8{x, r) = So sin[k ■ [x — e Ufca/a)], where 
e is the direction of the bulk velocity. Therefore 6 acquires a phase shift: 8cf){k) = k ■ e a/ a vj,. 
On squaring, and averaging over angles between k and e, this gives: 

{[6<P{lr)f) = l^,ev!{r). (53) 

Note that the averaging over angles is consistent with the angular integral done to get 
equation (52), and amounts to estimating the typical phase shift due to a superposition of 
bulk fiows of magnitude vj, directed randomly with respect to k. Substituting for v'^ from 
equation (13) and assuming it is dominated by the contribution from low A;, we recover the 
result in equation (52). Thus we have shown that for n < —1 the dominant phase shift 
is due to the kinematical effect of the bulk motion on small scales imparted by long wave 
modes. This is consistent with the connection between divergences in the nonlinear terms 
in the fiuid equations and vj, made in Section 2. 



4.2. Taylor Series Expansion 



In this section we make an additional assumption about the A;-dependence of 6[k) and 
6[k): we assume that in a small neighborhood around A;, 6 and 6 are smooth, differentiable 
functions of k. With this assumption we expand the nonlinear integrals in equations (8) in 
a Taylor series in (ki/k) about and restrict the range of integration to small ki. Thus 
we write: 6[k — ki) ~ 6[k) — ki ■ dS/dk^ and likewise for 6[k — ki). Unfortunately, the 
standard assumption about 6 and 6 in cosmology is that they are Gaussian random fields 
at the initial time. Thus at each value of k they are given by a random number drawn from 
a probability distribution. The distribution of 6 oi 6 with respect to k is quite the opposite 
of a smooth function, because its values at any two k are uncorrelated. We return to this 
point later in this section, but here we proceed with the Taylor series approach. 

With the Taylor expansion described above, the right-hand sides of equations (8), 
denoted by C{k) and D{k), take the form: 



C{k) = - / d-'k^ 



eCh) ^ Uk) - k ■ ^ j + e{k)8{k) + . . . 



(54a) 
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Dik) 



d^ki e{ki 



kl 



kjk^ 



dk 



(54b) 

where both equations have been expanded to the same order. In equations (54) we have 
included the contributions from [k — ki) ^ as well. We now write equations (54) at 
the order shown as linear equations by approximating 8 and 6 at small ki by the linear 
solutions 8i and 6i. Recall that we had obtained one linear term on the right-hand side 
of each equation in the previous section by introducing the integral a. Here we introduce 
three new integrals: ?y, 7 and (/ij, 



r, = d^ki6{ki) ■ 7 = d^ki[2{k ■ kf - l\8{ki) ■ = / d^k^kAjSik^) 



(55) 



where ku and kij denote the zth and jth components of the vector ki^ so that g^j IS a 
tensor. Note that aside from the dependence of 7 on the direction of A;, all the integrals in 
equation (55) are independent of k and r. In addition, all the integrals are convergent in 
an rms sense for n > — 3. 



We proceed by writing down a second order differential equation for 6 in terms of C 
and D: 

■■ n ■ fi .a 

(56) 



We then use the definitions of equation (55) to rewrite equations (54) as: 

C(A;) = iak ■ ab ^ a\kigi,jd,j\b — arjO , 

and 

D[k) = —iak ■ aO — d[kigijdj]6 — a^6 , 



(57) 



(58) 



where = d/dk^^ and the repeated indices i and j are summed over. We will now attempt 
to solve these equations for the amplitude and phase of (5 to a given order in a. We begin by 
using equations (57) and (58) to eliminate 6 in the terms on the right-hand side of equation 
(56) (we will also need to use the left-hand sides of equations (8)). Some algebra yields the 
following equation for 6: 



■■ a 

s + - 

a 



1 



1 -\- ar] 



7a — 2iak ■ a — 2a[kigijdj) 



6 

6 -(1 + ar])8 



a \k ■ a — i{kigijdj\n8)\ 6 — id \—kigijaj-\-i[kigijdj)ln6\6 



-i — • a — i{kigijdjln6) 



1 



1 -\- ar] 



7a + 



aa 



d? 



6 = 0. 



(59) 
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We now make a WKB analysis, which rehes on taking the phase to be more rapidly 
varying than the amplitude, to obtain self-consistent equations for the amplitude and phase. 



and g and F are 3 by 3 matrices, with / being the identity matrix. This solution can be 
verified by substitution using equation (59). Note that for ag < 1, F can be expanded as 
a Taylor series: F ~ a + a!^g l2 + a^g^ /6 + .... For ag <^ 1 the leading order solution is 
(f> = ak ■ a -\- (pi^ in agreement with (48). The solution for (f> in equation (61) can be used in 



Note that in this equation for A, all terms involving i and the divergent integral a have 
canceled out exactly! Hence the solution for A has no dependence on a, the only divergent 
integral in equation (59). Obtaining the full solution for A is still not possible as it requires 
solving equation (62), a nonlinear partial differential equation; however, for our purpose the 
key goal was the assessment of the a-dependence of A. Thus the Taylor series approach 
leads to two striking results: the solution for the phase given by equation (61), and the 
result that the evolution of the amplitude is not infiuenced by any divergent phase integrals. 

The above conclusions thus support the interpretation discussed in Section 2 that for 
—3 < n < — 1, there is no dynamically relevant divergence. Hence the evolution of scale free 
spectra will obey the standard self-similar scaling provided the statistics used are relevant 
to the growth of density perturbations. The divergent growth of the phase is not a measure 
of perturbation growth as it arises from bulk motions. However these conclusions rest on 
the assumption of a valid Taylor series approximation for 6 and 6. This assumption cannot 
be justified in the cosmological context for random-phase Gaussian initial conditions. It 
can be argued that the Taylor series expansion becomes a reasonable approximation at 
sufficiently late times when nonlinear evolution reduces the stochastic character of the 
density field. However this is at best a qualitative argument, and we must therefore regard 
the conclusions of this section as being suggestive of the answer, but not proven results. 



After some algebra we get the following relation for (f> by solving equation (59) at 0{a^): 




(60) 



This yields the solution. 
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In Paper II we shall measure the growth of the amplitude and phase of the density field 
in scale free simulations. The results from an n = —2 simulation in particular will allow us 
to test the analytical results of this section and to assess the validity of our approximations. 

5. Summary 

The goal of this paper was to examine the self-similar scaling of initially scale free 
cosmological spectra, P{k) oc k^. We emphasized that the scaling properties for n < —1 
have not been adequately studied either analytically or through N-body simulations. Indeed 
some results from N-body simulations suggested that the scaling properties of n = —2 
simulations were different from those with n > —1. In this paper we have used analytical 
techniques to investigate the possible breaking of self-similar scaling for n < —1. 

We motivated this investigation by demonstrating that terms diverging due to the 
contribution from long wave modes were present in the cosmological fiuid equations, but 
arguing that a detailed analysis was needed to assess their dynamical infiuence on the 
growth of density perturbations. In Section 3 we examined perturbative contributions to the 
power spectrum to examine the possibility of long wave divergences in these contributions 
for — 3 < n < — I. We found that divergent terms were indeed present, but that the leading 
order divergences canceled out at each order in perturbation theory. Terms which diverged 
less strongly were also present, but due to the computational complexity involved we did 
not proceed further with the perturbative analysis. We developed a non-perturbative 
approximation to study the nonlinear coupling of long wave modes in Section 4. We 
obtained a solution for the phase shift of the Fourier space density 6[k^ a) which is divergent 
for — 3 < n < — I. This divergence was interpreted as arising from the kinematical effect of 
the bulk fiows induced by long wave modes. With additional assumptions requiring that 
6[k^a) be amenable to a Taylor series expansion around A;, we showed that the evolution of 
the amplitude of 6[k^a) is not infiuenced by the divergent terms. It was emphasized that 
the Taylor series expansion cannot be justified for Gaussian random fields, and therefore 
the conclusion about the amplitude cannot be regarded as a proof. 

Thus the two approaches in this paper strongly suggest that the self-similar scaling of 
the amplitude of density perturbations holds for —3 < n < 4. The subtleties that arise 
due to long wave divergences for n < —1 affect statistical measures such as the rms phase 
shift, which does not measure dynamical nonlinearity. The general lesson is that certain 
statistical measures are susceptible to kinematical as well as dynamical infiuences, and must 
therefore be interpreted with caution. These include the rms particle displacement and any 
other statistic that involves the bulk velocity or measures driven by it. In Paper II we study 
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the scaling properties of scale free spectra in N-body simulations and make comparisons 
with the analytical results of this paper. 

It is a pleasure to thank Alan Guth for many stimulating discussions which helped 
to clarify and sharpen our arguments. We also acknowledge useful discussions with Neal 
Katz, Samir Mathur, Mark Metzger, Rajaram Nityananda, David Weinberg, Rien van de 
Weygaert and Simon White. 
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